/*
 * ibdscanner.cpp
 *
 *  Created on: Sep 27, 2011
 *      Author: lindenb
 */
#include <iostream>
#include <cerrno>
#include <cmath>
#include <cassert>
#include <fstream>
#include <cstdlib>
#include <set>
#include <vector>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <limits>
#include "xstdio.h"
#include "xstring.h"
#include "xzlib.h"
#include "throw.h"
using namespace std;

#define WHERE(a) cerr << __FILE__ <<":"<< __LINE__<<":" <<  a << endl


class Scanner
    {
    public:
	set<int> observed_columns;
	double max_diff;
	ostream* out;

	class Segment
	    {
	    public:
		std::string chrom;
			int start;
			int end;
			double p_value;
			int count;
	    };

	class Contigs
		{
		public:
			int col;
			string name;
			std::vector<Segment*> segments;
			Contigs()
				{

				}
			~Contigs()
				{
				while(!segments.empty())
					{
					delete segments.back();
					segments.pop_back();
					}
				}
			void print(std::ostream& out)
				{

				for(size_t i=0;i< segments.size();++i)
					{
					Segment* seg=segments[i];
					assert(seg!=NULL);
					if(i==0 || segments[i-1]->chrom.compare(seg->chrom)!=0)
						{

						seg->start=0;
						}

					if(i+1==segments.size() || segments[i+1]->chrom.compare(seg->chrom)!=0)
						{

						seg->end=numeric_limits<int>::max()-10;//-10 avoid side effect later
						}

					out << name << "\t"
						 << seg->chrom << "\t"
						  << seg->start <<"\t"
						  << seg->end <<"\t"
						  << seg->count <<"\t"
						  << seg->p_value
						  << endl
						  ;
					}
				}
		};

	vector<Contigs*> contigs;

	Scanner():out(&cout)
	    {
	    max_diff=0.1;
	    }

	~Scanner()
		{
		while(!contigs.empty())
			{
			delete contigs.back();
			contigs.pop_back();
			}
		}

	void run(Reader* reader)
	    {
	    vector<string> tokens;
	    vector<string> headers;
	    char* line;

	    if((line=reader->readLine())==NULL) THROW("Cannot read first line");
	    if(line[0]!='#') THROW("Expected first line to start with #");
	    XCString::split(line,"\t",headers);
	    for(set<int>::iterator r=observed_columns.begin();
	    		    r!=observed_columns.end();
	    		    ++r)
			{
			if((*r)>= (int)headers.size()) THROW("column="<<(*(r)+1)<<" out of range in "<< line);
			Contigs* ctg=new Contigs;
			ctg->col=(*r);
			ctg->name=headers[*r];
			contigs.push_back(ctg);
			}

	    *out << "#pair\tchrom\tstart\tend\tcount\tpvalue-obs" << endl;
	    while((line=reader->readLine())!=NULL)
			{
			if(line[0]==0) continue;

			XCString::split(line,"\t",tokens);
			if(tokens.size()!=headers.size()) THROW("Runtime error number of column!=header in "<< line);



			string curr_chrom(tokens[1]);
			int curr_position=atoi(tokens[2].c_str());



			size_t i=0;
			for(set<int>::iterator r=observed_columns.begin();
				r!=observed_columns.end();
				++r)
				{
				if((*r)>= (int)tokens.size()) THROW("column="<<(*(r)+1)<<" out of range in "<< line);
				Contigs* ctg=contigs.at(i);
				assert(ctg->col==(*r));
				char *p2;
				double obs=strtod(tokens[*r].c_str(),&p2);
				if(*p2!=0) THROW("bad number "<< *r << ":"<< tokens[*r]);
				if(	ctg->segments.empty() ||
					ctg->segments.back()->chrom.compare(curr_chrom)!=0 ||
					fabs(ctg->segments.back()->p_value - obs) > max_diff
					)
					{
					Segment* segment=new Segment;
					segment->chrom=curr_chrom;
					segment->start=curr_position;
					segment->end=curr_position;
					segment->p_value=obs;
					segment->count=1;
					ctg->segments.push_back(segment);
					}
				else
					{
					ctg->segments.back()->count++;
					ctg->segments.back()->end=curr_position;
					}
				++i;
				}
			}

		for(size_t i=0;i< contigs.size();++i)
			{
			contigs[i]->print(*out);
			}

	    }
    };

static void usage(std::ostream& out,int argc,char** argv)
		{
		out << argv[0]<< "\nLast compilation:"<< __FILE__ << " "<< __LINE__
		<<"\nPierre Lindenbaum PhD.\n"
		<<"Options:.\n"
		<<" -o (fileout)\n"
		<<" --obs (observed column index). This option can be used several times.\n"
		<<" --diff (val) Max difference between two p-values.\n"
		;
		out << endl;
		}

int main(int argc,char** argv)
    {
    Scanner scanner;
  
    int optind=1;

    char* fileout=NULL;
    fstream* fout=NULL;
    while(optind < argc)
	    {
	    if(std::strcmp(argv[optind],"-h")==0)
		{
		usage(cout,argc,argv);
		return EXIT_SUCCESS;
		}
	    else if(std::strcmp(argv[optind],"-o")==0 && optind+1<argc)
			{
			fileout=argv[++optind];
			}
	    else if(std::strcmp(argv[optind],"--obs")==0 && optind+1<argc)
			{
			int n=atoi(argv[++optind])-1;
			if(n<4)
				{
				THROW("Bad obs column \""<< argv[optind]<< "\"");
				}
			scanner.observed_columns.insert(n);
			}
	    else if(std::strcmp(argv[optind],"--diff")==0 && optind+1<argc)
			{
			scanner.max_diff=atof(argv[optind]);
			}
	    else if(std::strcmp(argv[optind],"--")==0)
			{
			++optind;
			break;
			}
	    else if(argv[optind][0]=='-')
			{
			std::cerr << "unknown option " << argv[optind] << std::endl;
			usage(cout,argc,argv);
			return (EXIT_FAILURE);
			}
	    else
			{
			break;
			}
	    ++optind;
	    }

   if(scanner.observed_columns.empty())
   	{
   	cerr << "Obs columns missing\n";
   	usage(cout,argc,argv);
   	return EXIT_FAILURE;
   	}

    if(fileout!=NULL)
	    {
	   fout=new fstream(fileout,ios::out);
	    if(!fout->is_open())
		    {
		    cerr << "Cannot open \""<< fileout << "\" " << strerror(errno) << endl;
		    return EXIT_FAILURE;
		    }
	    scanner.out=fout;
	    }





    try
	    {
	    if(optind==argc)
		{
		GZReader input;
		scanner.run(&input);
		}
	    else if(optind+1==argc)
		{
		GZReader input(argv[optind++]);
		scanner.run(&input);
		input.close();
		}
	    else
		{
		cerr << "Illegal number of arguments." << endl;
		return EXIT_FAILURE;
		}
	    }
    catch(std::exception& err)
	    {
	    cerr << "An error occured:" << err.what()<<  endl;
	    }
    catch(...)
	    {
	    cerr << "An error occured" << endl;
	    }
    if(fout!=NULL)
	    {
	    fout->flush();
	    fout->close();
	    scanner.out=&cout;
	    }
    return EXIT_SUCCESS;
    }
